/**
 *  DaisyMath
 *  Author: thanhkm
 *  Description: a model mathematic of Daisyworld
 */

model DaisyMath

/* Insert your model definition here */

global {
	
	// it is used for grow liminosity gradually
	int nb_time <-0;		
	
	float step_luminosity <- 0.001;
		
	//TEMPERATURE OPTIMAL for daisy to grow
	float TEMP_OPTIMAL <- 22.5;// = 295.5°K
	// FACTOR
	// Solar limunosity which changes in time
	float Solar_Luminosity <- 0.6;//+(time*(1.2/200));
	// Stefan-Boltzmann constant W/°K^4
	float SB_constant <- 5.67e-8;
	//W/m2 -- for reference, our Sun cranks out 1370 W/m2
	float Solar_Flux_Constant <- 864.65; 
	//this controls how the local temperatures of the daisies differ from the average planetary temperature
	float heat_absorp_fact <- 20.0;
	
	// ALBEDO
	// albedo of daisy and ground
	float black_albedo <- 0.25;
	float white_albedo <- 0.75;
	float ground_albedo <- 0.5;
	// Albedo of planet which will change on running
	float planetary_albedo <- 0.5;
		
	// Surface portion
	float Uncovered_Area <- 1.0; //portion of surface that daisies can occupy
	float Black_Area <- 0.0;
	float White_Area <- 0.0;
	
	// TEMPERATURE
	// Planet with daisy
	float Planet_Temp <- 253.0 - 273.0; 
	// Planet with no dasisies
	float Dead_Planet_Temp <- 253.0 -273.0; 
	// Temperature of black and white
	float Temp_Black_Land <- heat_absorp_fact*(planetary_albedo-black_albedo)+Planet_Temp;
	float Temp_White_Land <- heat_absorp_fact*(planetary_albedo-white_albedo)+Planet_Temp;

	// GROWTH RATE
	// the growth rate of Black and White
	float Black_Growth_fact <- 1-0.003265*((TEMP_OPTIMAL-Temp_Black_Land)^2);
	float White_Growth_fact <- 1-0.003265*((TEMP_OPTIMAL-Temp_White_Land)^2);		
	// death rate of black and white
	float death_rate <- 0.3;
	init{
		create DaisyMath number:1{
			
		}
	}
	reflex stop_simulation when: Solar_Luminosity > 4.0 {
		do halt ;
	}
}

entities {
	species DaisyMath{

		reflex growth{
			// Solar luminosity, it changes on the fonction of time: step_luminosity/cycle 
			Solar_Luminosity <- Solar_Luminosity + step_luminosity;
//			nb_time <- nb_time + 1;
		}
		
		reflex update{
			// Surface portion
			Uncovered_Area <- 1.0 - Black_Area - White_Area; //portion of surface that daisies can occupy
			
			// Albedo of planet
			planetary_albedo <- (Uncovered_Area*ground_albedo)+(Black_Area*black_albedo)+(White_Area*white_albedo);
					
			// TEMPERATURE
			// Planet with daisy
			Planet_Temp <- ((Solar_Luminosity*Solar_Flux_Constant*(1-planetary_albedo)/SB_constant)^0.25)-273; 
			// Planet with no dasisies
			Dead_Planet_Temp <- ((Solar_Luminosity*Solar_Flux_Constant*(1-0.5)/SB_constant)^0.25)-273; 
			// Temperature of black and white
			Temp_Black_Land <- heat_absorp_fact*(planetary_albedo-black_albedo)+Planet_Temp;
			Temp_White_Land <- heat_absorp_fact*(planetary_albedo-white_albedo)+Planet_Temp;

			// GROWTH RATE
			// the growth rate of Black and White
			Black_Growth_fact <- 1-0.003265*((TEMP_OPTIMAL-Temp_Black_Land)^2);
			White_Growth_fact <- 1-0.003265*((TEMP_OPTIMAL-Temp_White_Land)^2);		
			
		}

		// Daisy math equation
		float t;
		float h <- 0.1;
		equation DaisyEquation{
			// Equation of growth
			// the 0.001 is needed to give the system a bit of a nudge; without it the daisies never get going
			diff(Black_Area,t) =  Black_Area*(Uncovered_Area*Black_Growth_fact-death_rate)+0.01;
			diff(White_Area,t) =  White_Area*(Uncovered_Area*White_Growth_fact-death_rate)+0.01;
		}
	    solve DaisyEquation method: "rk4" step: h cycle_length: 1/h ;
		
	}
}

experiment DaisyMaths type: gui {
		parameter "Solar luminosity: " var: Solar_Luminosity category: "Planet" step: 0.1;
		parameter "Change of luminosity/cycle: " var: step_luminosity category: "Planet" step: 0.0001;
		parameter "Albedo of ds black: " var: black_albedo category: "Daisy" step: 0.02;
		parameter "Albedo of ds white: " var: white_albedo category: "Daisy" step: 0.02;
		parameter "Albedo of ground: " var: ground_albedo category: "Daisy" step: 0.02;

	output {
		monitor "luminosity" value: Solar_Luminosity; 
		display Daisyworld_charts refresh_every: 1 {
			chart "World temperature evolution" type: series background: rgb("white") size: {1, 0.4} position: {0, 0.05}{
				data "world temperature without daisy" value: Dead_Planet_Temp color: rgb("blue");
				data "world temperature with daisy" value: Planet_Temp  color: rgb("red");
				}
				
			chart "Daisy areas cover" type: series background: rgb("white") size: {1, 0.4} position: {0, 0.55}{
				data 'Black area' value: Black_Area color: rgb('black') ;
				data 'White area' value: White_Area color: rgb('blue') ;
				data 'Ground area' value: Uncovered_Area color: rgb(216, 74, 6) ;
			}
		}
		
		display Solar_luminosity refresh_every: 1 {
			chart "Solar luminosity" type: series background: rgb("white") size: {1, 0.4} position: {0, 0.25}{
				data 'Solar luminosity' value: Solar_Luminosity color: rgb("blue") ;
			}
		}
	}
}
